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Reduced density matrix hybrid approach: An efficient and accurate method 
for adiabatic and non-adiabatic quantum dynamics 
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Department of Chemistry, Stanford University, 333 Campus Drive, Stanford, California 94305, 
USA 

We present a new approach to calculate real-time quantum dynamics in complex systems. The formalism is 
based on the partitioning of a system's environment into "core" and "reservoir" modes, with the former to be 
treated quantum mechanically and the latter classically. The presented method only requires the calculation 
of the system's reduced density matrix averaged over the quantum core degrees of freedom which is then 
coupled to a classically evolved reservoir to treat the remaining modes. We demonstrate our approach by 
applying it to the spin-boson problem using the noninteracting blip approximation to treat the system and 
core, and Ehrenfest dynamics to treat the reservoir. The resulting hybrid methodology is accurate for both 
fast and slow baths, since it naturally reduces to its composite methods in their respective regimes of validity. 
In addition, our combined method is shown to yield good results in intermediate regimes where neither 
approximation alone is accurate and to perform equally well for both strong and weak system-bath coupling. 
Our approach therefore provides an accurate and efficient methodology for calculating quantum dynamics in 
complex systems. 
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I. INTRODUCTION 

The calculation of real-time quantum dynamics for 
large molecular systems is a longstanding problem in 
chemical physics. Of particular interest is the case of 
a small subsystem embedded in a surrounding thermal 
bath, which forms the basis for the investigation of con- 
densed phase energy and electron transfer as well as spin 
and charge transport in nanoscale devicesjii^ In such sit- 
uations, one is typically interested in the calculation of 
the system's reduced dynamics, averaged over the bath 
degrees of freedom. Due to the importance of these 
problems and the absence of a general solution, a vari- 
ety of methods have been developed which vary as to the 
regimes in which they are accurate and their applicability 
to large systems. Many of these approaches rely on either 
averaging over semi-classical trajectories of bath degrees 
of freedom or propagating the system's reduced density 
matrix (RDM) directly without explicit treatment of the 
bath, for example by master equations or path integral 
techniques. 

Trajectory based methods include the mean-field 
Ehrenfest, '^'4 surface-hopping^i^ semi-classical initial 
value representationjIiS. linearization)2r22 and classical 
mappingi^"— approaches. These offer the appeal of 
a transparent link to the classical limit and simple 
application to complex condensed phase systems and 
to slow baths. However, these methods typically fail 
in regimes with strong system-bath coupling or high- 
frequency baths where quantum effects are known to be 
important. 

Techniques based on direct propagation of the re- 
duced density matrix include the Markovian and non- 
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Markovian Bloch-Redfield equations^®"^^ and the nonin- 
teracting blip approximation^'^^ which serves as a time- 
dependent generalization of Forster theory To make 
such approaches manageable, perturbation theory in the 
system-bath coupling (Bloch-Redfield) and in the non- 
adiabatic coupling (NIBA, Forster) are typically carried 
out only to second-order and hence performance degrades 
when these terms become large. In addition, path inte- 
gral approaches, such as quantum Monte Carlc^^"^'^ and 
the quasi-adiabatic propag ator path integral (QUAPI>2i 
have also been developed to propagate the RDM. While 
these approaches can be made numerically exact, they 
are frequently difficult to converge in practice when the 
system-bath coupling is strong or the bath is slow. Fur- 
thermore, many variants scale poorly with the number 
of system states, drastically limiting the size of treatable 
systems. 

Hence when the bath dynamics are "fast" compared to 
those of the system (the so-called non-adiabatic regime), 
both perturbative quantum master equations and numer- 
ically exact path integral approaches are typically accu- 
rate since the correlation time of the bath is short, render- 
ing non-Mar kovian memory effects less significant. When 
the bath dynamics are "slow" compared to those of the 
system (the adiabatic regime), these methods fail but it 
becomes reasonable to treat the bath degrees of freedom 
classically and hence trajectory-based schemes generally 
yield accurate results. Given this complementarity, it 
would clearly be desirable to develop a scheme which can 
accurately treat both regimes by naturally tuning be- 
tween the above approaches in a consistent manner. In 
addition, such an approach might effectively treat inter- 
mediate regimes where no obvious time-scale separation 
exists and hence neither RDM nor trajectory-based ap- 
proaches are accurate. 

The "dynamical hybrid" scheme of Wang, Thoss and 
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Miller—!^ offers such an approach in which bath de- 
grees of freedom are partitioned into "core" modes to 
be treated with quantum mechanics and "reservoir" 
modes to be treated with classical mechanics. In this 
scheme, the wavefunction-based multi-configurational 
time-dependent Hartree (MCTDH)^"— method is used 
for the system and core while a mean-field Ehrenfest 
treatment is applied to the reservoir. When applied to 
the spin-boson problem, this dynamical hybrid scheme 
frequently converges to the exact result with only a small 
percentage of modes included in the core. However, to 
achieve convergence in the presence of fast baths, up to 
80% of the bath modes need to be included in the core. 
Due to the high cost of treating modes using MCTDH, 
this severely reduces the efficiency of the dynamical hy- 
brid approach and hence limits the ease with which it 
can be applied to larger systems. 

In this work, we present a methodology which avoids 
the use of wavef unctions, replacing them with a sin- 
gle reduced density matrix averaged over the quantum 
core modes. The reservoir modes exert a classical time- 
dependent driving force on this system RDM, which in 
turn drives the classical modes via a back-reaction force. 
This reduced density matrix hybrid formalism (RDM- 
Hybrid) allows one to use an appropriately chosen par- 
titioning of the bath into core and reservoir modes such 
that the system RDM can be propagated either by an 
exact reduced dynamics algorithm with a decreased cost 
or by an inexpensive approximate master equation with 
improved accuracy. In other words, by using physical 
intuition, one can choose the core-reservoir partitioning 
such that different methods are applied only to the set of 
bath modes for which they are expected to work. 

We demonstrate our approach by applying it to the 
spin-boson model and treating the system and core us- 
ing the noninteracting blip approximation (NIBA}iiiS 
This approach yields a methodology which is numerically 
cheap while obtaining impressive quantitative accuracy 
over all regimes investigated as well as offering excellent 
scaling with the number of system states. 

The outline of the paper is as follows. Section III Al 
briefly reviews the Ehrenfest methodology. This back- 
ground is then used in Sec. IIIBl to derive the RDM- 
Hybrid approach. Section IIIII then describes how this 
scheme can be applied to the spin-boson model. The 
results of the application are presented in Sec. IIVI and 
Sec. IVl concludes. 



II. THEORY 

A. Ehrenfest dynamics 

While there are many routes to the mean-field Ehren- 
fest equations of motion and subsequently the RDM- 
Hybrid method to be presented below, we proceed via 
the quantum-classical Liouville equation, which provides 
a particularly simple derivation. 



We begin by considering the generic coupled system- 
bath Hamiltonian, 



H = H,{s,ps) + i/fc(Q,P) + Vsbis, Q) 
where the free bath Hamiltonian is 



and the coupling is assumed to be of the form 
K6 = 5]j^.(s)G,(Q). 



(1) 
(2) 

(3) 



The dynamics of the total density matrix, p(t), is given 
by the Liouville equation. 



dpjt) 
dt 



i[H,p{t)], 



(4) 



where henceforth we set h ^ I. By employing the partial 
Wigner transform of the density matrix with respect to 
the / environmental degrees of freedom. 
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(27r)/ 



d^ 



ye 



X (Q + ||/5|Q-|), (5) 



and taking the classical limit of the corresponding partial 
Wigner transformed Liouville equation (for more details, 
see e.g. Ref. Isol ). one obtains the so-called quantum- 
classical Liouville equation 



+ '-{{H^,p^it)} + {p^it),H^}) 



(6) 



where [•, •] denotes a commutator and {■, ■} is the Poisson 
bracket. 

Now, one makes the approximation that the system 
and bath reduced density matrices (RDMs) decouple at 
all times, 



p^(Q,P,i) = Ps(t)Pfc(Q,P,t), 



(7) 



with the system RDM ps{t) = J d^Qj dfPp^{Q,P,t) 
and bath RDM pb(Q,P,t) = Tr^ys [p^(i)] . Inserting 
this product form into the quantum-classical Liouville 
equation and noting that the bath reduced density matrix 
is simply the classical phase space distribution. 



pt{Q,F,t)^S{Q~Q{t))S{P-Pit)) 



(8) 



one obtains the mean-field Ehrenfest equations of 
motionSi for the system RDM, 



dt 



= -t[Hs + Vsb{s,Qit)),Ps{t)] 

= -z[i7,+^i^,(s)G.(Q(t)),Ps(t)] 

i 

= -i[Hs{t),ps{t)], 



(9) 
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and the classical bath degrees of freedom, 
dQk 



dt 
dPk 
dt 



Pk 



(10) 



d 



V{Qk)+'Tr,{Vsb{s,Cl{t))] 



dV{Qk) ^dG,{Q{t)) 



dQk 



dQk 



TIs{F^Ps{t)}. (11) 



The Ehrenfest approximation to an observable of the 
reduced system is given by 

(0,(t)) = Tr,Trb{03p(i)} « Tr^jTr, {O,p,(0}pt(t)} 
~ Jd^^J dfPTTs{0,ps{t)}pbiQ,P,t). (12) 

In accordance with the classical treatment of the bath, 
this latter integral is evaluated by molecular dynamics 
obeying the equations of motion given above. To gener- 
ate the canonical ensemble average consistent with the 
initial density matrix 



p(0) = p.(0)p6(0) = psiO)exp{~l3Hb)/Zb, 



(13) 



one averages over classical trajectories with bath initial 
conditions sampled from either the Wigner-transformed 
bath density operator. 



pr(Qo,Po) = 



1 



{2n)f 



df 



ye 



X (Qo + ||/5fc|Qo-|), (14) 



or more simply from the classical Boltzmann distribution. 



Pb 



■(Qo,Po) = e~''^^('^-^°V^. 



(15) 



Although the Wigner distribution gives initial condi- 
tions consistent with the exact quantum distribution, this 
property is not conserved by the classical dynamics. 



B. Reduced density matrix hybrid dynamics 

In contrast to the approach taken above, in which the 
full Hamiltonian is split into a system and bath with 
coupling between the two, we instead divide the bath 
degrees of freedom into two sets: the reservoir modes 
{{Qk,Pk ■ k = 1,...,/'}) to be treated classically and 
the core modes {{qk,Pk ■ k — f' + l, ■ ■ ■ , /}) to be treated 
quantum mechanically. Hence, we partition the Hamil- 
tonian as 

H = Hsc{s,Ps, ci, p) + HriQ, P) + Vsr{s, Q) (16) 

where the system-core Hamiltonian is 
Hsc{s,Ps,<i,p) = Hs{s,ps) + Hc{q,p) + Vsc{s,q) (17) 



and the system-core and system-reservoir couplings are 
defined as 



Kc(s,q)=^i^.(s)G,(q), 

i 

14r(s,Q) -^^^.(s)G,(Q). 



(18) 
(19) 



Since nothing in our above discussion of Ehrenfest dy- 
namics utilized the specific properties of the system, we 
could just as well re-label the combined system and core 
as the system of the previous section. Likewise, we as- 
sociate the reservoir with the bath from before. Thus 
we assume the total (Wigner-transformed) density ma- 
trix factorizes at all times into a system-core RDM and 
reservoir RDM, 



p'^{t)=Psc{t)pb{Q,P,t). 



(20) 



Taking the classical limit we arrive at the Ehrenfest-like 
RDM-Hybrid equations of motion for the system-core 
density matrix. 



dpscjt) 
dt 



-l[Hsc + V,r{s,Clit)),Psc{t)] 

= -I [Hs + J2 (Q(^)) ' 

i 

= -^[7^,(^),p,,(^)], (21 

and the classical bath degrees of freedom, 
dQk 



dt 

dPk d 

dt 



Pk (22) 
V{Qk) + TTsTr,{V,b{sM{t))} 



dQk 
dV{Qk) 
dQk 



2^ J^. Tr,Tr4J^,p,,(t)|. 



dQk 



(23) 



Although the exact calculation of the system-core den- 
sity matrix dynamics given by Eq. (|2ip is extremely 
demanding, this was essentially the approach taken by 
Wang, Thoss, and Miller (WTM)25.'26 averaged high- 
dimensional wavefunction trajectories with initial condi- 
tions of the core wavefunction sampled from the exact 
quantum mechanical core density matrix and classical 
reservoir degrees of freedom sampled from the quantum- 
classical Wigner distribution. 

The important point we seek to make in this work 
is that Fi is a pure system operator, i.e. the core and 
reservoir are not directly coupled, such that the classical 
reservoir equation of motion above may be written 



dPk 
dt 



dVjQk 
dQk 



where we have introduced the system RDM, ps{t) = 
Trc{psc(^)}, averaged over the quantum core degrees of 
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freedom. Furthermore, the calculation of any dynamical 
system variable is given by an expression analogous to 
the pure Ehrenfest approximation above, 

{0,{t)) = Tr.TreTr, {O.pit)} « Tr,{Tr, {OsPs{t)} Pr{t)] 
« J df'Q J df'PTT,{0,p,it)}pr{Q,P,t), (25) 

which again only requires the system RDM (and not the 
much higher-dimensional system-core RDM), with clas- 
sical reservoir trajectories to be sampled from an appro- 
priate distribution. 

We have thus arrived at a self-consistent dynamical 
scheme which treats the core bath modes quantum me- 
chanically yet only requires the propagation of a typical 
system RDM and classical reservoir degrees of freedom. 
This RDM-Hybrid formulation presented above consti- 
tutes the major theoretical result of this work. We point 
out that a similar approach appears to have been inves- 
tigated by Golosov et al.^^ though it was restricted to 
their previously developed memory-equation algorithnt^^ 
and only minimally pursued. 

When the system RDM averaged over the quantum 
core modes, Ps{t), is treated fully quantum mechanically, 
our methodology yields the exact result when all bath 
modes are included in the core. However, in practice 
one expects to achieve numerical convergence before this 
limit is reached. Of course, the exact quantum dynamical 
treatment of the system RDM is in general no less nu- 
merically challenging than the original problem (a system 
coupled to a quantum bath). However, the advantage of 
this approach is that one can tailor the partitioning in 
Eq. (fT6|) either to alleviate the computational expenses of 
a numerically exact reduced quantum dynamics method 
(such as the QUAPI-based iterative tensor propagation 
scheme of Makri et al.^^-'^^"'^^) or to improve the accuracy 
of approximated quantum dynamics (such as a pertur- 
bative quantum master equation) for the system RDM. 
In this sense, the above presented formalism outlines a 
framework for the efficient partitioning and calculation of 
real-time reduced quantum dynamics where one can em- 
ploy any exact or approximate dynamical scheme for the 
system RDM depending on the complexity of the prob- 
lem. As discussed in the following section, here we adopt 
an approximate treatment of the system RDM using a 
perturbative master equation to yield a very cheap hy- 
brid methodology with impressive quantitative accuracy. 



III. APPLICATION TO THE SPIN-BOSON PROBLEM 

While the above formalism is indeed quite general here 
we demonstrate its effectiveness at treating the spin- 
boson Hamiltonian, 

H = ea,+Aa,+a, ^ CkQk+Y, [Pk + ^IQI] /2, (26) 

fe k 



which describes a two-level system with energy bias 2£ 
and tunneling matrix element A, linearly coupled to the 
coordinates of a bath of harmonic oscillators. The spin- 
boson Hamiltonian has been used as a model for a variety 
of physically distinct quantum relaxation and transport 
processes. 

A principal observable of interest is the population 
variable, 

P{t) = TrsysTrbath Wz{t)pm = {a,{t)), (27) 

which measures the difference in population of sites 1 
and 2. Furthermore, we will henceforth work in reduced 
units with respect to the tunneling matrix element. A, in 
addition setting h = ks — ^- To begin our application, 
we first discuss the separation of modes into the core and 
the reservoir. 



A. Separation of modes 

The harmonic bath in the spin boson problem is fully 
specified by its spectral density 

k ^'^ 

When separating the bath modes into a core and reser- 
voir, we will require that 

Jcorc(w) + Jros(w) J(w). (29) 

Intending to treat high frequency modes in the quantum 
mechanical core and low frequency modes in the classical 
reservoir, we use a switching function, S{ll!,uj*), which 
switches from 1 to as a; goes from to w*. This allows 
the core and reservoir spectral densities to be defined as, 

Jcorc(^) = J(^)[1-5(C^,W*)], (30) 
J,es(t^) = J(t^)5(w,W*). (31) 

An infinite number of switching functions satisfy such a 
criteria and the best form will depend on the Hamiltonian 
adopted, the methodology used to treat the core, and 
the required balance been accuracy and efficiency. The 
simplest approach would be to use a step function, 

Siu.,.n = \l (32) 
10 uj > uj . 

However in many cases, such a sharp switch will leave 
a core which is expensive to treat using quantum me- 
chanics. The reason for this difficulty can be seen by 
considering the bath correlation function, given by 

k 

1 

— — dujj{uj)[coth{(3uj/2)cos{ujt) — iam{ujt)] , (33) 
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CO* = 4Q3 



(o* = 8(0 




FIG. 1. The effect of splitting tiie spectral density (top) on the 
resulting bath correlation functions with T/luc = 4 (bottom). 



whose correlation time determines the range of memory 
(non-Markovian) effects. It is straightforward to see that 
using a step function switching hke that in Eq. p2|) will 
result in a long-time oscillatory tail in the bath corre- 
lation function, a(t), and thus extensive non-Markovian 
effects which are difficult to treat quantum mechanically. 
It is therefore advantageous to separate the spectral den- 
sity using a smooth switching function which can yield a 
short bath correlation time for the core modes and allow 
for efficient quantum mechanical treatment. 

For our purposes a switching function of the form, 



B. Noninteracting blip approximation for the system and 
core 

Although the RDM-Hybrid framework allows for the 
application of many approaches to treat the system and 
core modes, in this work we employ the noninteracting 
blip approximation (NIBA)iiii^ To see why such a choice 
is computationally simple, note that because the fluc- 
tuating classical degrees of freedom provide an effective 
time-dependent bias, the system Hamiltonian no longer 
commutes with itself at different times, and thus its prop- 
agator is a cumbersome time-ordered exponential (Dyson 
series). However, since NIBA only treats the diagonal 
part of this Hamiltonian exactly, its propagator is triv- 
ially calculated. 

Another important reason for our choice is the comple- 
mentarity between NIBA and Ehrenfest. While Ehren- 
fest dynamics work best in the adiabatic regime, Wc/A <^ 
1, NIBA is perturbative in A/ujc, thus working best in 
the non-adiabatic regime, ujc/A. ^ 1. To exemplify this 
complementarity, in Fig. [2] we show a coherent to in- 
coherent crossover, as a function of Wc/A in the high- 
temperature, strong-coupling regime. As can be seen in 
Fig. [21 the Ehrenfest method yields quantitatively exact 
population dynamics for an adiabatic choice of parame- 
ters, Wc/A <^ 1 [panels (a) and (b)], whereas NIBA gives 
the exact result upon crossing over into the non-adiabatic 
regime, Wc/A ^ 1 [panels (c) and (d)]. With this in mind 
we are strongly encouraged to consider a NIBA-Ehrenfest 
hybridization, the details of which are described next. 



C. Simulation details 



The total bath was discretized by the relation. 



S{lu,lu*) = 



1 - iLu/uj*Y 



UJ < LU 



LO > OJ 



(34) 



was chosen. In Fig.[Tl we show the bath correlation func- 
tion when this function is applied to the Debye spectral 
density. 



J(w) = 2XuJc 



UJ 



(35) 



where A — n^^ dojJ{uj)/uj is the reorganization en- 
ergy. As can be seen, this smooth switching does not 
introduce any long-time tails and allows the bath corre- 
lation time for the core modes to be decreased as w* is 
increased and more of the slow modes, which give rise to 
temporally non-local effects, are moved into the reservoir. 
In addition, removal of slow modes from the core encour- 
ages the use of approximate master equations, such as 
the noninteracting blip approximation discussed in the 
next section, which complement the Ehrenfest treatment 
of the slow modes. 
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2 J{uk) 

-^k—, T: 



(36) 



where p{uj) is a density of frequencies chosen to reproduce 
the reorganization energy for any number of modes, /. 
For an Ohmic spectral density J(w) = Xuj / uJcF{lo /ujc) 
with cutoff function F{lj/uJc), we thus require 



TT 



duj 



k=l 



A F{ujk/ujc) 



p(Wfc) 



which is clearly satisfied by the choice 



/ 



TTUJc 



(37) 



(38) 



For the Debye spectral density considered here^ the cutoff 
function is given by F{u!/ujc) = 2/[l + (uj/uJc) ]■ 
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FIG. 2. Population dynamics, P{t) = {(^zit)), for the spin- 
boson model with a Debye spectral density. The Hamiltonian 
parameters are e = 0, A = 2.5A, and T = 2A. Ehrenfest 
and NIBA are compared to the numerically exact results of 
Thoss, Wang, and MillerSS in (a)-(c) and to our own con- 
verged QUAPI calculation in (d). 



Having obtained a discrete distribution of frequencies, 
we then take into account the switching function and 
recalculate the reservoir couplings as 



2 Jrcsi^k) 
Ck = -Wfc— T — ^, 
TT p{UJk) 



(39) 



in the process removing all bath modes with loj. > w* 
from the reservoir. In a complimentary fashion, the core 
spectral density, Jcorc(w), is used in the NIBA equations. 

In all results presented, / = 300 discrete bath modes 
were found to be sufficient for convergence and averages 
were performed over 10^ — 10* reservoir initial conditions. 
For pure Ehrenfest results, we sample from the Wigner 
distribution, which for the spin-boson Hamiltonian's har- 
monic bath takes the form. 



/ 



pr(Q,p) = n 



tanh(^Wfc/2) 



k=l 



X exp 



2tanh(/3c^fc/2) f ^ u;lQl 



(40) 



As alluded to above, when high frequency baths are 
present, Wigner sampling can improve the accuracy of 
Ehrenfest dynamics since the initial conditions are con- 
sistent with the exact quantum mechanical distribution. 
However, since the zero-point energy inserted by Wigner 
sampling is not conserved by the classical dynamics, long- 
time predictions can be unreliable (zero-point energy 



leakage) Hence, in regimes with strong coupling or 
high bath frequencies, sampling from either the classical 
Boltzmann or quantum Wigner distributions can produce 
inaccurate results. In contrast, for our RDM-Hybrid ap- 
proach we find that the results are largely insensitive to 
the choice of reservoir initial conditions, since most quan- 
tum mechanical modes are included in the core leading 
to only low frequencies being present in the reservoir, 
for which classical and Wigner sampling give identical 
results. This effect is shown later, in Fig. 21 where pure 
Ehrenfest dynamics using classical sampling is seen to un- 
derestimate the decay of the system population variable 
while using Wigner sampling overestimates the decay. 
The RDM-Hybrid approach, however, produces graph- 
ically identical results using either sampling scheme. 

The time evolution of the coupled classical dynamics 
and NIBA equations were solved with a second-order 
Runge-Kutta scheme with a timestep of O.OIA"^. For 
the spin-boson Hamiltonian considered here, the reser- 
voir equations of motion, Eqs. take the form 



dQk 
dt 

dPk 
dt 



Pk 



= -i^lQk 



CkP{t), 



(41) 
(42) 



with the system population variable P{t) — 
Tr^ {cr2/9s(i)}. Taking P{t) to be constant over a 
half-timestep as called for by the Runge-Kutta scheme 
used here, one has 



Pk{t + At/2) 



Qu{t) + ^P{t) 
^k 

.^^sm{ujkAt/2) 



cos(wfeAt/2) 



''\p{t), (43) 



Qk{t) + ^P{t) 



Pk{t) 



cos(wfeAi/2) 



Wfe sin(a;feAt/2) 
(44) 



For a two-level system coupled to a common bath, the 
NIBA master equation for the population difference with 
a time-dependent bias is given simply as^ 



dP{t) 
dt jQ 

with the kernels 



dT 



K+{t,T)P{T)+K^{t,T) 



(45) 



K+{t, t) = 4A2 exp [-Q2{t - r)] cos [Qi (t ~ r)] 

cos{C(t,T) + {1 + S) [Qi(t) - gi(T)]} , (46) 

K_ {t, t) = 4A2 exp [-Q2{t - r)] sin [Qi(< - t)] 

sin{C(i, r) + (1 + 6) [Qi(t) - Qi(r)]} . (47) 

The parameter S originates from the bath initial con- 
dition, with (5 = if the system and bath are initially 
uncoupled and S = —1 if the bath is initially solvating 
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the donor. The twice-integrated bath correlation func- 
tion, Q{t) = Q2{t) + iQi{t), is given by 



4 r°° 

Q{t) = - doj 
Jo 



■h^I^ {coth(;9w/2) [1 - cos(wi)] 



-|- i sin(wt)} . 



(48) 



The effect of a general time-dependent bias 2£{t) is 
captured in the function 



C(i,T) = 2 / dt'e(t'), 



(49) 



which for the core-reservoir splitting considered here is 
defined as 



dt' 



ckQk{t) 



k G res 



(50) 



The true population difference, P{t) = {<Tz{t)) is finally 
given by the average of P{t) over classical trajectories 
with initial conditions sampled from the appropriate dis- 
tribution. 

Importantly, the scheme described here yields the cor- 
rect population dynamics of a two-level system driven 
by a classical external field, so that by taking the limit 
Lu* — >■ cxD, we remove all bath modes from the core, and 
recover the Ehrenfest solution. In the opposite limit, 
a;* — 0, we recover the full NIBA treatment. 

Because we are not treating the core exactly, it is im- 
portant to emphasize that uj* is no longer a convergence 
parameter. Instead, u* tunes the balance between the 
two approximate methods, here taken to be Ehrenfest 
dynamics and NIBA. The optimal a priori choice of this 
switching frequency is an interesting problem for future 
work, but here we present a physically motivated pre- 
scription for its determination. We again recall that 
"adiabaticity" is inherently a problem of timescale sepa- 
ration. The adiabatic regime is physically realized when 
the timescale of the bath is much greater than that of 
the system, and vice versa for the non-adiabatic regime. 
A given bath mode k could be classified into "core" or 
"reservoir" by comparing its characteristic timescale, oj^T^ 
with that of the system, u!~yg. For the spin-boson Hamil- 
tonian, Eq. (l26t . the system frequency Wgys is given by 
the Rabi frequency. 



A2. 



(51) 



In practice we use the smooth switching function, 
Eq. (|34|). and following the line of argument above, set 
the switching frequency equal to the Rabi frequency, i.e. 
a;* — ujR. This procedure effectively partitions the bath 
modes into those that are faster than the system dynam- 
ics (comprising the core) and those that are slower (com- 
prising the reservoir). While comparison with existing 
exact results has shown that this choice, w* — ujr, is not 
always truly optimal, we will show in the next section 
that it nevertheless provides a very robust methodology, 
with quantitative predictive ability. 



(a) e = 0, = 5A, co = 0.25A, T = 2A 



I » Exact 

RDM-Hybrid 

Ehrenfest 

NIBA 




2 4 6 8 10 12 14 16 
(b) e = 0, X = 2.5A, co = 0.25A, T= 0.2A 




10 12 



Time t [A 



FIG. 3. Population dynamics of the spin-boson model in the 
adiabatic regime (a;c/A < 1) without an energetic bias. The 
bath has a Debye spectral density and parameters as given in 
the figure. Note the difference in the time axis. 



IV. RESULTS 

In all our results we compare to the numerically ex- 
act population dynamics computed by Thoss, Wang, and 
Miller^^ for the spin-boson Hamiltonian, Eq. (pS)) . with 
a Debye spectral density, Eq. (j35p and initial condition 
p(0) = |1)(1| exp{-l3Hb)/Zb, i.e. (5 = in Eqs. (|46l)-(|4ll). 

We begin in the adiabatic regime, oJc/^ < 1, where 
Ehrenfest dynamics are known to be relatively accurate. 
Fig. [HJa) shows that this is indeed the case for high tem- 
perature, T = 2A. Though qualitatively good, the very 
strong coupling, A — 5A, degrades the accuracy of the 
Ehrenfest approach. However, our RDM-Hybrid method 
yields a long-time population decay in much better agree- 
ment with exact results. In Fig. [3l^b), we investigate an- 
other system with strong coupling, A = 2.5A, but with a 
temperature one order of magnitude lower. The classical 
mechanical approximation made by using Ehrenfest dy- 
namics intrinsically worsens for lower temperature, where 
quantum mechanical effects are expected to play a more 
dominant role. We see that the RDM-Hybrid result is 
again in better quantitative agreement with the exact 
long-time limit. Both panels (a) and (b) in Fig. [3] un- 
derscore the important point that our RDM-Hybrid ap- 



8 



e = 0, ?L = 5A, CO = 5A, T = 2A 

C 




o RDM-Hybrid (Classical IC) 
Ehrenfest (Wigner IC) 



-0.2 - ° Ehrenfest (Classical IC) 
NIBA 

-0 4' ■ ' ' ' ■ ' ' ' ■ ' ' ' 

2 4 6 8 10 12 

Time f [A"^] 

FIG. 4. Population dynamics of the spin-boson model in the 
non-adiabatic regime {ujc/ A > 1) without an energetic bias. 
The bath has a Debye spectral density and parameters as 
given in the figure. Also shown are the effects of classical vs. 
Wigner sampling of the bath (reservoir) initial conditions for 
Ehrenfest dynamics (RDM-Hybrid). 

proach is not simply an "average" of the two methods, 
but can in fact generate qualitatively different dynamics 
which are not situated "in between" those of the individ- 
ual methods. 

The system trapping seen above leading to a quasi- 
stationary state is nearly impossible to capture with per- 
turbative quantum master equation approaches, as ex- 
emplified by the poor NIBA results. For comparison, the 
WTM MCTDH-Ehrenfest scheme required 25% of the 
bath-modes to be treated quantum mechanically, sup- 
porting our observation of largely classical bath dynam- 
ics. 

In Fig. 21 we investigate the non-adiabatic regime, 
ujc/A > 1, with strong coupling A = 5A, where NIBA 
is quantitatively accurate. Here we show the effects of 
sampling from a classical distribution or a Wigner distri- 
bution. The Ehrenfest approximation with Wigner initial 
conditions yields a decay which is too rapid, while that 
with classical initial conditions yields a decay which is 
too slow. Both completely miss the two-step relaxation 
dynamics at very short times. RDM-Hybrid captures 
these short time dynamics exactly and shows a robust 
insensitivity to the initial sampling employed. Although 
our hybrid method exhibits a slight discrepancy at longer 
times, it is much more accurate than the Ehrenfest ap- 
proach. 

Having demonstrated that our RDM-Hybrid approach 
can yield very good agreement with the exact results 
when only one of its composite methods (Ehrenfest 
dynamics or NIBA) is successful, we now consider in 
Fig. [5] the intermediate regime, Wc/A — 1, where nei- 
ther method is particularly accurate by itself. Fig. ^a.) 
depicts the population dynamics at the high temperature 
T ~ 2A for relatively strong coupling, A — 2.5A. Here, 
both NIBA and Ehrenfest dynamics are in very poor 
numerical agreement with the exact result, though the 



(a) e = 0, = 2.5A, co^ = A, T = 2A 




2 4 6 8 10 12 14 16 18 
Time t [A^] 



FIG. 5. Population dynamics of the spin-boson model in the 
intermediate regime {uc/A = 1) without an energetic bias. 
The bath has a Debye spectral density and parameters as 
given in the figure. Note the difference in the time axis. 

Ehrenfest dynamics correctly predicts a two-step relax- 
ation, albeit only qualitatively. On the contrary, RDM- 
Hybrid correctly exhibits a rebound in the population 
dynamics near t ~ 2A-1 (although one that is somewhat 
overpronounced) and yields dynamics at all later times 
in excellent agreement with the exact result. 

In the lower panel, Fig.[5l^b), we consider a weaker cou- 
pling and a temperature two orders of magnitude lower. 
Here, NIBA does quite poorly and Ehrenfest dynamics 
are mildly better, exhibiting slightly overdamped oscilla- 
tions and a weak phase-shift at long times. Once again, 
RDM- Hybrid is nearly exact, despite this severely low 
temperature. 

In both examples above, i.e. Figs. [SJa) and (b), 
Thoss, Wang, and Miller found it necessary to treat 
over 50% of the bath modes quantum mechanically with 
MCTDH. It is remarkable that our computationally in- 
expensive RDM-Hybrid method employing NIBA for the 
core modes is able to yield such similarly accurate re- 
sults for these problems, which serves as a testament to 
its robustness. 

We conclude this section by briefly considering the ef- 
fect of an energetic bias on the performance of our RDM- 
Hybrid method. In Figs.[ni[a) and (b), we show two sets 
of population dynamics for the same high-temperature. 
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(a) e = 0, >L = 0.25A, o) = 0.25A, T = 2A 



(a) e = 0, ?L = 0.25A, co = 5A, 7= 2A 
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(b) e = A, X = 0.25A, oo = 0.25A, 7= 2A 
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FIG. 6. Population dynamics of the spin-boson model in the 
adiabatic regime (tUc/A < 1) without an energetic bias (a) 
and with an energetic bias (b). The bath has a Debye spec- 
tral density and parameters as given in the figure. Note the 
difference in the time axis. 
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{h)^ = A,X = 0.25A, co = 5A, 7 = 2A 




Time t [A" ] 



FIG. 7. Population dynamics of the spin-boson model in the 
non-adiabatic regime (wc/A > 1) without an energetic bias 
(a) and with an energetic bias (b). The bath has a Debye 
spectral density and parameters as given in the figure. 



adiabatic regime, but with an energetic bias e = A in 
panel (b) . As expected, NIBA performs poorly due to the 
adiabaticity of the conditions under consideration in both 
panels. As for the Ehrenfest method, while it gives excel- 
lent results in the unbiased case [panel (a)], it strongly 
deviates from the correct dynamics in the presence of 
a bias [panel (b)]. This deficiency of the Ehrenfest ap- 
proach and similar quantum-classical methods for biased 
system is well-known and typically connected to unphysi- 
cal treatment of zero-point energy in classical treatments. 
Attempts to correct this behavior have included an ad- 
justable zero-point energy parameter for sampling the 
initial conditions'^^ and its determination by comparison 
with exact short-time moments 4^ 

In spite of this potential concern, we see in Fig.[6]that 
the RDM-Hybrid approach performs incredibly well and 
equally so for both unbiased and biased systems. Unfor- 
tunately, this success is not universal, as seen in Fig. [71 
which considers the effect of a bias on non-adiabatic dy- 
namics. While RDM-Hybrid is very successful for the un- 
biased case, it performs more poorly in the presence of a 
bias. This error appears to be attributable to the Ehren- 
fest dynamics, which RDM-Hybrid tracks quite closely 
at short times. Nonetheless, the RDM-Hybrid popula- 



tion does appear to approach the correct long-time value, 
which both Ehrenfest and NIBA incorrectly predict. 

V. CONCLUSIONS AND FUTURE WORK 

To summarize, we have formulated a quantum-classical 
methodology for the quantum dynamics of a coupled 
system-bath Hamiltonian, averaged over the bath degrees 
of freedom, by employing a core-reservoir partitioning of 
the bath degrees of freedom. Unlike the technique of 
Wang, Thoss, and Miller;^ our RDM-Hybrid approach 
avoids the use of expensive wavefunction-based quantum 
mechanics, replacing it by a reduced density matrix cal- 
culation. This allows for flexibility to use a variety of 
inexpensive approximate approaches for different parts 
of the computation, allowing for combined accuracy and 
efficiency. 

Within the above framework, we then used the nonin- 
teracting blip approximation (NIBA) for the reduced dy- 
namics, yielding an efficient, scalable quantum dynamics 
methodology with excellent accuracy for both adiabatic 
and non-adiabatic parameter regimes. Specifically, for 
Ng discrete states, a numerical implementation of NIBA 
scales only as N"^ (see Eq. (25])), which allows for very effi- 



10 



cient investigation of bigger systems. The largest discrep- 
ancy with exact resuhs was found in the non-adiabatic 
regime with an energetic bias. To alleviate such deficien- 
cies, we enumerate a number of avenues for future work. 

One could imagine pursuing alternative approximate 
treatments of the quantum core-averaged system RDM. 
For example, the accuracy of NIBA (employed here) is 
known to degrade at low temperature, in particular for 
biased systems. The Redfield equations, on the other 
hand, have been very successfully applied to such sys- 
tems, being restricted, however, to weak-coupling and 
non-adiabatic baths. Thus, an RDM-Hybrid approach 
employing Redfield and Ehrenfest dynamics may be one 
way to extend the validity of the Redfield equations into 
these regimes where it would otherwise fail. 

Ideally, one would like an exact treatment of the quan- 
tum system RDM, such that the method becomes numer- 
ically exact when all modes are included in the core. As 
alluded to at the end of Sec. Ill B[ we believe the tensor 
propagation scheme of Makri et ali^i^"— is especially 
promising. Because the QUAPI algorithm scales expo- 
nentially with the memory time required to span the bath 
correlation function, the partitioning discussed here lead- 
ing to a reduced bath correlation time would constitute 
an enormous reduction in computational expense. Work 
along this direction is currently in progress. 

Lastly, for a chosen system-core dynamics method, a 
more rigorous investigation of the form of the switching 
function, S{uj,uj*), and switching frequency, w*, should 
be pursued. Ideally these choices should be automated 
and not left up to physical considerations as was done 
here. For example, one might try comparing the short- 
time moments of the population dynamics to those ob- 
tained exactly from quantum perturbation theory. The 
investigation of this point will be a subject of future work. 

In spite of these desired improvements, the RDM- 
Hybrid method employing NIBA has been shown to be 
very accurate, often correcting the long-time populations 
of the otherwise accurate Ehrenfest dynamics. Future 
work will include the application and investigation of 
both of these methods to models of electronic energy 
transfer in molecular aggregates, including the multi- 
site Frenkel exciton Hamiltonian describing the Fenna- 
Matthews-Olsen complex, a prototypical photosynthetic 
system. 
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